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Abstract 

We consider time-harmonic linear elasticity equations in domains 
containing two-dimensional semi-infinite strips. Since for such 
problems there exist modes with different signs of group and phase 
velocity, standard perfectly matched layer (PML) as well as standard 
Hardy space infinite element methods fail. 

We apply a recently developed infinite element method for a physi¬ 
cally correct discretization of such waveguide problems which is based 
on a Laplace transform in propagation direction. In the Laplace do¬ 
main the space of transformed solutions can be separated into a sum 
of a space of incoming and a space of outgoing functions where both 
function spaces are certain Hardy spaces. The Hardy space is chosen 
such that the construction of a simple infinite element is possible. 

The method does not use a modal separation and works on in¬ 
tervals of frequencies. On those intervals the involved operators are 
frequency independent and hence lead to linear eigenvalue problems 
when computing resonances. Numerical experiments containing con¬ 
vergence tests and resonance problems are included. 
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1 Introduction 


Computational methods for wave equations bear great attention due to their 
huge importance for real live problems. In solid mechanics the application 
reaches from simulating seismic waves, non-destructive testing to material 
characterization [10, 19]. Such equations are usually posed on unbounded do¬ 
mains fi. Mesh based methods, i.e. finite difference/volume/element meth¬ 
ods deal with that difficulty by truncation of hi to a bounded subdomain fl int 
and special treatment of the exterior domain f2 ext = \ 0; n t- 

Based on representation formulas for the solution in h2 e xt with given 
boundary data at the artificial boundary T = fi ext fl Hint local approxima¬ 
tions of the Dirichlet-to-Neumann operator at T can be used (see EH for a 
review). Also based on representation formulas for the solution (via a Green 
function) boundary element methods typically lead to non-local but more 
accurate approximations. 

A method, which does not directly use a representation formula, is the 
complex scaling method, reintroduced by Berenger [6] as perfectly matched 
layer method (PML) for electromagnetic waves. The method became soon 
very popular and can be classified as todays standard method for treating 
unbounded domains. In in a Hardy space infinite element method (HSIE) 
was introduced for Helmholtz problems. The method can be understood 
as a special infinite element method, relying on the pole condition [16] . It 
features many advantages of PMLs, in particular it does not depend directly 
on a representation formula. However, the theoretical background of this 
method is different to that of the PML. 

As for all linear wave equations, a fundamental technique in understand¬ 
ing linear elastic wave equations is the study of most simple solutions EM- 
In general unbounded domains these are plane waves e wt ' x_la,t w, where as 
in waveguides, such as plates R 2 x /,/ C R and cylinders Rxf),Dc R 2 , 
boundary conditions have to be respected by the solutions. The general form 
in these cases are modal waves. In semi-infinite cylinders they take the form 
e ira_lw# w(y), where t G R>o is the time variable, x G R>o is the longitudinal 
coordinate, y G D the vector of transverse coordinates, oj G R>o the angular 
frequency, k the wavenumber and w/|w| the direction of displacement. The 
frequency c o and the wavenumber k have to fulfill a wave equation specific 
dispersion relation oj(k) to yield a solution. The wave travels in direction 
sign(/-s;) with phase velocity oj/k, whereas the energy is transported in direc¬ 
tion sign(<9u;(/c)/chv) with group velocity du o{n)/dn [ 21] . 

The possibility of waves with different signs of group and phase veloc¬ 
ity was already discussed by Lamb in 1904 [T%] , Such mismatches between 
the directions of propagation of phase and energy cannot only happen in 
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waveguides, but also in anisotropic materials |5]. Although PML methods 
select waves by their propagation direction of phase, the attention to mis¬ 
matches between propagation directions of phase and energy in elastic ma¬ 
terials was omitted for long [5]. This may result in a wrong selection of 
outgoing/incoming waves and can lead to an exponential growth of solutions 
in the damping layer of PMLs. Such cases are usually reported as insta¬ 
bilities, due to the definition of stability, but first of all, in such cases an 
unphysical radiation condition is incorporated in the PML formulation. 

Sometimes, such problems can be solved using a transformation of vari¬ 
ables jH HJ. In [29, [22] SMART layers were introduced for seismic waves. 
They allow more flexibility in the damping than PMLs, but lose the prop¬ 
erty of being perfectly matched. For semi-infinite linear elastic cylinders in 
the time-harmonic setting, two methods based on bi-orthogonal relations of 
modal solutions were proposed: The method presented in [3] is based directly 
on a modal representation of the solution, whereas in [28, 8] ways to modify 
standard PMLs are reported. They rely on a smart post processing by ex¬ 
changing backward incoming with backward outgoing modes. Since, these 
methods involve modal solutions, which depend non-linearly on c o 2 , they 
depend themselves non-linearly on u A Thus if they are used to compute res¬ 
onances, they lead to non-linear eigenvalue problems. Although solvers for 
non-linear eigenvalue problems exist (e.g. 0 ). they are more involved than 
solvers for linear eigenvalue problems. 

In [13] a new family of Hardy space infinite elements was introduced. It 
is based on a generalized pole condition, which allows for different signs of 
phase velocities. The achievement of this paper is to pick up these results 
and apply it to two-dimensional time-harmonic semi-infinite cylindrical linear 
clastic waveguide problems. A big advantage of this method is, that u> 2 only 
enters as a scalar coefficient in the method. Hence, the discretization of a 
resonance problem leads to a generalized linear matrix eigenvalue problem, 
which can be treated with a standard solver. Moreover, the numerical 
results indicate a super-algebraic convergence with respect to the number of 
degrees of freedom in the longitudinal direction. 

The outline of the paper is as follows: Section [2] formulates the diffrac¬ 
tion problem to solve, in particular the modal radiation condition and its 
reformulation as a pole condition. In Section [3] we introduce the Hardy space 
infinite element in one dimension as in CHI and discuss the choice of method 
parameters for the investigated elasticity problem. Section [4] explains the 
discretization of the elasticity problem with the use of tensor product ba¬ 
sis functions. Section [5] deals with the spectral objects and properties of 
resonance problems. Finally, in Section [6] we give numerical examples. 
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Figure 1: sketch of an elastic waveguide problem under consideration in this 
paper 


2 General setting 

2.1 Geometry and elasticity Equations 

We are considering in this paper elastic waveguide problems in two dimen¬ 
sions (see Fig. [I] for a typical situation): The domain of interest 12 consists of 
a bounded interior domain 12 int and L semi-infinite waveguides W \,..., Wl 
with interfaces Ti,..., T^. The two-dimensional time-harmonic isotropic 
linear elasticity problem is given by 

— div <t(u) — poj 2 u = f in 12, (2.1a) 

B u = g on <912, (2.1b) 

u — u inc satisfies a radiation condition in each waveguide ll-y, i — 1,..., L. 

(2.1c) 

Here, 9ft(u(x)e~ lt ^) for x G 12 and time t > 0 is the time-harmonic displace¬ 
ment vector, p > 0 the density, co > 0 the angular frequency, A, p > 0 
are Lame parameters, e(u) = )(Vu + (Vu) T ) the strain tensor, cr(u) = 

A div u • Id + 2yue(u) the stress tensor, and f a volumetric force. In this paper 
we will always use the plain strain model (see e.g. 0), i-e- the Lame param¬ 
eters are related to Young’s modulus E and Poisson’s ratio u by p — yyyyy 
and A = Boundary conditions at <912 int D <912 are formulated in 
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terms of a trace operator B , e.g. a Dirichlet trace operator Bn = u or a 
Neumann trace operator Bn = cr(u) • n with outer normal vector n, and a 
boundary datum g. 

We will assume in the following, that p, A and p, are constant in each 
waveguide Wg , that there exists no volumetric force f in the waveguides, 
and that we have traction free boundary conditions er(u) • n = 0 at the 
boundaries dWg Ddfl. The terms radiation condition and incident field u mc 
in (2.1c) and in Fig. [l] will become clear in the next subsection. 


In addition to the scattering problem (2.1), where the angular frequency 
u and the sources f, g, and u mc are given, we will consider the corresponding 
resonance problem: Find resonances w 6 C with positive real part and non¬ 
trivial resonance functions u such that 

— div cr(u) — pu 2 u = 0 in 12, (2.2a) 

Hu = 0 on <912, (2.2b) 

u satisfies a radiation condition in each waveguide Wg, i = 1, 

(2.2c) 


,,L. 


2.2 Radiation condition 


In order to define a physically correct radiation condition (2.1c), we consider 
in this section one single waveguide. The generalization to multiple waveg¬ 
uides is straightforward. Let W := M + x (—R,R) be a reference waveguide 
and p, A and p be constant and positive. As for most other problems the 
radiation condition can be derived by an analytic representation of solutions. 
For waveguides a convenient representation is a modal sum. For given angu¬ 
lar frequency w>Owe call u(», «;u;) a mode with wavenumber k(oj) G C, if 
it has the form 


u(£,p;w) = e lK( ^ c w(p;o;), (£,p) G W,u G M + , (2.3) 

and solves 

— div <r(u) — pu 2 n = 0, (£, rj) G W, (2.4a) 

o"(u) • (?) = 0, (£, V) G M+ x {-R, R}. (2.4b) 

It is straightforward to see that if ^ ™) * s a m °de, so is ^ 

as well as e~ lK ^^ ) an d e XK ( u ^ ^ ) • A physical solution should 

be bounded for £ —>■ oo. Hence, if A(«:(a;)) ^ 0 we want to exclude modes 
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with A(k(o;)) < 0. In this case, we call the exponentially decaying modes 
evanescent. If A(k(o;)) = 0, both modes with wavenumbers ±k(lu) stay 
bounded and it is not obvious which one is physically relevant. We give here 
two approaches leading to the same distinguishing criterion. For both we 
assume, that 3?(9 w «;(u;)) flz 0, 9f(9 w /c(u;)) = 0. 


1. 121 states that the velocity of energy transport of a mode is given by the 
group velocity d(jj(n)/dK which has the same sign as dhi(uj)/duj. Since 
in problem (2.4) there is no source in W, energy should be radiated to 
infinity and therefore dn(oj)/du> > 0. See also [17., Rem. 3.1]. 


2. For the limiting absorption principle we add an artificial small damping 
e > 0 to the system, i.e. we substitute lu > 0 by lu + ie. The corre¬ 
sponding solution u e to the damped version of ( 2.4[ ) should be bounded 

8k(u) £ 


for x —> 00 . Since by linearization | exp (in(u> + ie)£) | ~ exp 
this leads to dn((ju)/d(ju > 0. For further details see PT7| Cor. 3.1]. 


duj 


Therefore, if dn{pj)/duj > 0 we call e lK ^ ( ) an outward propagating 


mode and e ^ j an inward propagating mode. Note, that there 

exist modes with different signs of phase and group velocity (see Fig. [2j) . We 
call a wavenumber k(lu) outgoing , if there exists an evanescent or outward 
propagating mode of the form (2.3) satisfying (2.4). In other words k(cu) is 
outgoing (incoming), if 

1. d u n(uj) > 0 (d u K(uu) < 0) for real wavenumbers k(l 0 ) E M, and 


2. A(k(c<;)) > 0 ( 9 (^( 0 ;)) < 0) for non-real wavenumbers k(co) $ M. 

A function satisfies the radiation condition in the waveguide W and is called 
outgoing , if it can be approximated by a linear combination of evanescent 
and outward propagating modes. It is not a priori clear, that the traces 
of evanescent and outward propagating modes on the waveguide interface 
T = {0} x (— R, R ) are dense in L 2 (Y) 2 . In the following Remark we cite 
some results from m, which justify our definition of the radiation condition. 


Remark 2.1 In [77j / a rigorously analysis of modal decompositions for two 
and three dimensional elastic waveguides is given using a quadratic eigen¬ 
value problem. There, not only modes (eigenfunctions) of the form (2.3) are 
considered, but also the generalized eigenspaces of the investigated quadratic 
eigenvalue problem (see H3 (0.10)]). We refer to these functions of the gen¬ 
eralized eigenspace, which are not eigenfunctions, as associated modes in this 
paper. We briefly summarize the main results of 
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Figure 2: left: the first eight dispersion curves. Modes corresponding to the 
red solid part have positive group velocity and modes corresponding to the 
blue dashed part have negative group velocity, right: the first nine outgo¬ 
ing/incoming wavenumbers multiplied with i marked with red squares/blue 
diamonds at frequency cu = 1.615. Parameters are H = p = E = 1, and 

i/ = 0.2. 


!■ E3 Thru. 1.5]: For fixed to > 0 there exist wavenumbers K n (uj),n G N, 
symmetrically situated relative to the real axis and the origin. They are 
situated in arbitrarily small angles, adjoining the imaginary axis, with 
the exception of a finite number of wavenumbers. In particular, only a 
finite number n n (uj) G M. exists. 

2. 23 Thm. 3.8]: There exists a sequence of frequencies 0 = uj\ < oj\ < 

■ ■ ■ < co 2 —)■ oo, such that for all to G M + \ {co n , nGN} the implication 
K n (ui) GR=t- d u K n (uj) Gl \ {0} holds. 

3 . 23 Thm. 2.5]: Assume u j G M + \ {co n ,n G N}. Then the traces of 
modes and associated modes corresponding to {n n {ijj), n G FT: ^(hi n (co)) > 
0 V (k„(w) G RA d u hi n (co) >0)} are dense and minimal in L 2 ( Y) 2 as 
well as in H 1 ( Y) 2 . 


2.3 Numerical implementations of the radiation con¬ 
dition 


In a numerical method for (2.1) the radiation condition as to be taken into 
account. Due to the simple waveguide geometry, methods based on modal de¬ 
compositions as the one proposed in [3] are popular. Assuming that every so- 
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lution in a waveguide can be expanded in a sum of modes ^ ngN e lKn ^w n (r)), 
a Dirichlet-to-Neumann operator can be defined as 


DtN^i 




nGN 


W «( ? ?) : = y^CTn 

nGN 


( e iK "^w, 


(v))- 


(2.5) 


Truncating the waveguides W[ and posing cr n (u) = DtN u at the artificial 
ends leads to a formulation of (2.1) in a bounded domain. For a numer¬ 


ical method standard finite element methods for the bounded domain can 
be used together with a discrete Dirichlet-to-Neumann operator (DtN)#, 


which can be constructed by using only a finite number of modes in (2.5), 
i.e. (DtN)# E„<jv eiK " HCw «W := lt can be shown 

that the modes satisfy biorthogonal relations, which can be exploited to im¬ 
plement (DtN)# in an elegant way, as done in [3]. However, the method 
requires a precomputation of the modes. Moreover, since the modes depend 
on u, the discretization matrix will depend non-linearly on oj. Hence, dis¬ 


cretizing the resonance problem (2.2) leads to a large, non-linear eigenvalue 
problem. 

In m a numerical method for a Helmholtz problem is proposed based on 


eigenfunction expansions of the bounded interior problem coupled to mode 
expansions in the waveguides. These approaches lead again to non-linear, 
but comparatively small eigenvalue problems. Nevertheless, solving the non¬ 
linear eigenvalue problem numerically is a non-trivial task and of course a 
computation of the modes and interior eigenfunctions is needed. 

A popular method which does not depend on the computation of the 
modes is the perfectly matched layer method. The longitudinal direction £ 
of the waveguide is thereby complex scaled by a G C, i.e. £ := a£, leading to 
a new equation for u(£, rj) := u(£, rj). If a is chosen such that 3?(ai/c n (o;)) < 0 
for all outgoing wavenumbers K n (u) and u is an outgoing solution of the orig¬ 
inal problem, than u is exponentially decaying in £. Hence, a truncation of 
the infinite waveguide to a finite waveguide with zero Neumann or Dirich- 
let boundary condition at the end introduces only a small error due to the 
exponential decay. The resulting equation can then be discretized with stan¬ 
dard methods, as it is posed on a finite domain. The linear system matrix 
takes the form K ^ — aAM# When looking for resonances, the above leads to 
generalized linear matrix eigenvalue problems, which is a main advantage in 
view of numerical algorithms. 

The reason why perfectly matched layers are only of limited use for waveg¬ 
uide problems with backward propagating modes, is that no a G C can exist, 
such that SR(ak n (w)) < 0 for all outgoing wavenumbers K n (u) (see Fig. |2(b) , 
a multiplication with a leads to a rotation of i K n in the complex plane). 
The modified perfectly matched layer methods presented in [28], 8j combine 








complex scaling with a special treatment of the backward propagating mode. 
For scattering problems with only a few backward propagating modes this 
approach works well, but again it leads to non-linear eigenvalue problems 
when discretizing the resonance problem (2.2). 


2.4 Pole Condition 

Our goal is to construct a radiation condition, which combines the advan¬ 
tages of the methods presented in the previous subsection: It should yield 
physically correct solutions and discretizations of resonance problems should 
lead to linear matrix eigenvalue problems. Again for the sake of simplicity, 
we formulate the following only for the reference waveguide W. For co > 0 
let 

N 

u(f, r l) = '52 Cn U n (f, V), (f, V) e W, 

n 

with c n G C be a finite sum of modes and associated modes with outgoing 
wavenumbers K n (oj) G C. Then u is bounded and the Laplace transform in 
the longitudinal direction £(u(», rj))(s) for ( s,r /) G M + x (—R,R) takes the 
form 


N 


£(u(«,77))(s) = 


V Kr, 




( 2 . 6 ) 


with m n G N. It has a meromorphic extension to C with poles at {ift n }. In 
contrast, the Laplace transform of e~ lKn % has a pole at — \K n (see Fig. 2(b) for a 
typical situation). The pole condition states, that u is called outgoing, if £u 
has no poles in a suitable region of the complex plane. Here, one might chose 
a point symmetric, smooth and asymptotically straight boundary curve T 
separating the outgoing from the incoming poles, simply connected domains 
T~, r + such that C = r - urur + and asks £ u to be holomorphic in T - . 
Typical choices of F can be seen in Fig. [4j The assumption 


c r H 


(2.7) 


with §(cj) := (i/ci(o;), i/C 2 (w),... } being the set of outgoing wavenumbers 
multiplied with i for the frequency u, is thereby essential. To formulate 
the above in a more compact way, we introduce the Hardy space H~(T ) as 
the subspace of L 2 (r), such that for each / G H~{T) there exists an in T - 
holomorphic function / vo i with / being the non-tangential limit of f vo We 
refer to m App. A] for a detailed definition and properties of such Hardy 
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spaces. We can formulate the pole condition now in the following way 


£u e [r(r) ® l 2 (T)] 2 . 


( 2 . 8 ) 


If (2.7) holds true, then the modal radiation condition of Sec. 2.2 


is 

equivalent to the pole condition ( |2.8[ ) for solutions u to Equ. ( |2.4[ ) having the 
form u = ^ 2^=1 c n u n and u n being modes and associated modes. 

The pole condition (2.8) does not depend on the wavenumbers but only 
on T and is therefore frequency independent. The wavenumbers are hid¬ 
den in the assumption ( 2.7[ ). Of course, we have to discuss how to chose 
T (see Sec. 3.3). Moreover, for the interpretation of solutions to the reso¬ 
nance problem (2.2) in Sec. [5] we need the wavenumbers. But the numerical 
method based on (2.8) presented in the following sections is independent of 
the wavenumbers and waveguide modes. This facilitates the implementation 
of the method a lot. Moreover, a discretization of the resonance problem 
(2.2) leads to a linear generalized matrix eigenvalue problem. 


3 The Hardy space infinite element 

The pole condition £u 6 [i7 _ (T) ® L 2 (Y)] 2 is a nice way to reformulate 
the physical radiation condition, but a stable numerical method based on 
this framework is delicate. In ms such a method was developed for a one 
dimensional toy problem including a complete convergence analysis. We 
present here only the results needed for an implementation in the setting of 
a convected one dimensional Helmholtz equation 

— u" + v! — u 2 u = 0, x > 0, (3.1a) 

m'(0) = u' 0 , (3.1b) 

CueH~(T), (3.1c) 


whereas the curve T is chosen such that the solution to (3.1) is unique. Since 
solutions of (3.1a) have the form C\e' K1 + Cie tK1 with C±, C2, «i, K2 G C, the 
curve T is such that ifiq G T + and i/v 2 € T~ (or vice-versa). 

The discretization of the elasticity problem discussed in Section [4] will 
be straightforward. In order to enhance readability, we have refrained from 
giving the correct mathematical framework including the variational formu¬ 
lation in the Hardy space. For the clastic waveguide problem this can be 
deduced along the lines of ra and ag¬ 


io 
























3.1 Basis functions and infinite element matrices 


Similar to a classical infinite element method we start with the variational 
form 

/»oo /»oo /»oo 

/ u'vdx'+ / u r vdx — uj 2 / uv dx =—u o v(0) (3.2) 


of (3.1), which holds true for the solution u and sufficiently fast decaying, 


smooth test functions v. We are using a Galerkin scheme with the same basis 
functions <^’ ong , j € N, as ansatz and test functions. These basis functions 
should fulfill the pole condition £<^ ong <E H~(T). 

From [133 we know that for any two complex parameters So,Si € T + the 
linear hull of 




SOiSl, 


s : = 


Sq+Si / S + Sq 


Si 


So 


Lb+i)/2j 


S + Si 
S — Si 


U/2J 


j E No, (3.3) 


is dense in the Hardy space H (T). The basis functions <^ ong of our Galerkin 
scheme are defined via their Laplace transforms 


(£^l ong )(s) := 


(r long\ / v ^-^(S) 

(£<Pi )(») : = 


s - So 


3 = 2, • 


s - So 
They have the form 

tp l ° ng {x) = e S0X pf(x ) + e slx p s Mx), x > 0, 


s e r. 
(3.4) 


with polynomials pf ,pf and in particular span{<^Q° ng , yy ong } = span{e So3; , e SlX }. 
Comparing the basis functions with the modes ( 2.3| ) shows, that these func¬ 
tions represent the longitudinal part of the modes exactly for wavenumbers 
ft = s 0 /i and ft = s\/i. Therefore, s 0 /i and si/i can be considered as two dif¬ 
ferent wavenumbers. The actual choice of these parameters will be discussed 
in the next subsections. 

It is easy to show that <^ ons , j E N fulfill the pole condition and moreover, 
by a limit theorem of the Laplace transform there holds 


ion g (o) _ ii m (^ long (x) = lim s(£^ ong )(s) = 

J :r\<0 J s—>oo J 


3 = 1 
3 > 1 


This facilitates coupling of these basis functions at x = 0 with e.g. standard 
finite element basis functions for x < 0. 

We are left with the computation of the integrals in (3.2) using y/° ng for 


u and v. If sq and si are chosen with negative real part, the integrals are 
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bounded. Similar to m Lemma A.l] it can be shown, that 

= g^(£^"«)( s )(£rf”')(- s ), ds, j,k<= N. 

(3.6) 

Since £y^ ons are meromorphic functions, the integrals on the right hand side 
can be computed by the residue theorem. Following the computations in 
[13] . the infinite mass matrix (M^) := / 0 °° </9’ OIls (x)</9[? ns (a;) dx is given by 

the tridiagonal matrix 


MZ g = 


SqSi 


Si -Si 


0 0 

— Si SQ+Sl 


-so 0 




0 

-SO 

0 

0 


V 


So + Sl 

-Si 

-Si 

So+Sl 


(3.6) 


Analogue calculations can be performed for d x ip^, since 

„ sib S0 ’ Sl (s) 

(Cd xV 'D(s) = (Cd,vY*)(a) = -im, 3 = 2,.... 

s- S 0 J S- So 

The drift (D™ ng ) jk := ,/|f d x <p l ° De (x)<p l ™ e (x) dx and stiffness (S£j jk := 
/ 0 °° 9 x <p 1 ° ng (x)d x ip 1 ^ ns (x) dx matrices are given by 


/ i 


i 


D 


OO 

long 


1 

2 


-i 


V 


o 

-i 


i 

0 


\ 



so so 

So So+Sl 


0 0 

si 0 

o o 

So+Sl So 

so SQ+Sl 


V 


(3.7) 

Using only the first Ad ong basis functions and therefore only the fV long x A long 
upper left block of these matrices leads to the discretization of (3.2): Find 
f/( Arl ° nR ) G C^° ng such that 


/ „jv lon g T-jlV lon s 
l ^long ' -^long 



(-u'o, 0 ,..., 0 ) 


T 


(3.8) 


Note, that the matrices 5' 1 ^ l | " r ' e , and M^° ns are independent of the 

angular frequency oj. 


3.2 Relation between the method parameters sq? si and 

r 


is a conforming discretization of (3.la)-(3.lb). The side constraint 
G H~(T) (see (3.1c)) is however not only fulfilled for 


- N l°n g (TV lo”*) long 

L - 2^ 7 = 1 U j Yj 
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the chosen T, but for any T with so, si G T + . If for such a second T the 


wavenumbers ft 12 of (3.1a) multiplied with i belong to different sides of the 


curves, e.g. iK\ G r + fir and i k 2 G F 
with the condition Cu G H 


flf + , Equations (3.la)-(3.lb) together 
(r) define a different solution then (3.1). But 


(3.8) is again a conforming discretization for this second problem. This leads 


(p l ° ng of (|3.8h converge 


m 


to the question that if the solutions Y2f=i 8 Uj N E ' ) ' 
any norm) to a function U, does C U G H~(T) ov CU G H~(T) hold? 

The answer was given in P! : The linear hull of the functions £y^ ong , j G 
N, is dense in any H~(T) with s 0) s i G F + . But they form a stable basis only 
in one specific space H~ (T soso ), i.e. any U G H~(T SOtS1 ) can be expanded U = 
a i ^ Fj° ng w ifh a square summable sequence (aq) je N- One can deduce 
that if {ifiq,bv 2 } O T ± = {i/Vi, i/v 2 } 0 Eq. ( 3.8| ) is uniquely solvable at 
least for sufficiently large N long with solutions and lim^viong^.^ || Cu — 


ST Nlong rr(iV lolls ) r .long 

Z^i= l u i ’-'Y 


3 ~ rj 11 (r so SJ ) = 0 holds. Thence the choice of s 0 ,Si for 
implicitly poses the pole condition w.r.t. T S0 S1 on the solution u of 


The curve T 


S0)Sl is given by the algebraic variety 

S + Sq s + Si 


r , 0)S1 < s G C: 


S — Sq S — 


= l 


(3.9) 


In order to have a physically correct pole condition for elastic waveguide 
problems we are left with the task to find s o, Si such that §(o;) C T+ Si , which 
will be the issue of the next subsection. The following properties of the curve 

|s-£o| |s-Sj| 


T SQiSi can all be deduced from (3.9): First of all for g SQiSl (s 


have the characterizations 


|s+so| |s+Sl| 


we 


r )Li = { s G C: 9s 0 , Sl (s) < !}, F S0)S1 = {s G C: g 80tai (s) = 1}, and T = { 


The curve can explicitly be parameterized by T SOiSl = 7 S0i 


7s 0 ,si( r ) := 


r 2 (s 0 + si) + |s 0 | 2 si + 

Si 

| 2 S 0 

k 2 (s 0 + Si) + |s 0 | 2 Si + 

Isll 

| 2 s 0 | 




r G 


with 


(3.10) 


(3.11) 


which shows that its asymptotic behavior is a line. In the case |s’o = |si| it 
actually is a line. Further scaling the poles scales the curve: r rs „. rjI1 = rT s 
for all r > 0. If we chose so, si such that 


3?(so),9ft(si) < 0, 
9(so + si) > 0, 
| So I 2 Q : (si ) + |si| 2 9(s 0 ) < 0, 


(3.12a) 

(3.12b) 

(3.12c) 
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^ C • 9sq,si (s) ' 



























CO 


Figure 3: The first eight dispersion curves for H = p = E = 1 and v = 0.25. 
Modes corresponding to the red solid part have positive group velocity and 
modes corresponding to the blue dashed part have negative group velocity. 


it holds i(—C(s 0 ,Si),0)Ui(C(so,Si),oo) C T+ and i(0, COo, Si))Ui(-oo, -C(s 0 , «i)) C 


r so , sl with 


COo,si) 


po| 2 Q(si)+|si| 2 3 ; (5o) 

Otso+si) 


(3.13) 


3.3 Choice of sq and s\ 


Coming back to the elastic waveguide problem, we inspect the location and 
qualitative properties of the wavenumber spectrum in order to find reason¬ 
able values for the parameters s o and si. First, let us look only on the 
real wavenumbers in Fig. [3] for fixed parameters v and H,p,E. There exist 
eight frequencies oj > 0 which support waves with vanishing group velocities 
(ihf) • them yield vanishing wavenumbers K n (uj) = 0 and two don’t. 

These 8 frequencies cui,... ,cu 8 are part of the sequence of frequencies defined 


in Rem. 2.1 and cannot be treated by the pole condition, since for cu —> ujj for 


j = 1 ,... ,8 an outgoing wavenumber converges to an incoming wavenum¬ 
ber. Hence, outgoing wavenumbers cannot be separated by T from incoming 
wavenumbers. In Fig. [3] we have made a (non unique) decomposition into 
five intervals (a, b ) such that for all c o G (a, b ) \ {a^!,... ,a; 8 } one of the two 
following cases hold true: 
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1. For all real wavenumbers the signs of group velocity and phase velocity 
coincide. Hence, the real wavenumbers with positive group velocity 
belong to (0, oo). 


2. There exists 9 > 0, such that the real wavenumbers with positive group 
velocity belong to (— 9, 0) U (9, oo). 

Hence, the first case requires i(0, oo) c T + and the second case i(— 9, 0) U 
i(9,oo) C r+. Chosing different parameters H,p,E would only result in a 
different scaling of Fig. [3j On the other hand varying v changes the qualita¬ 
tive behavior. Nevertheless, for all oo either the first case or the second case 
seems to be true for cylindrical elastic waveguide problems with constant E , 
p, and v in the waveguides (at least the authors were not able to find any 
contradicting parameters). 


For the non-real outgoing wavenumbers it is known (see Rem 2.1), that 
for a fixed frequency they are contained in a circular proper subsector of 
(s 6 C: 9s > 0}. Because the wavenumbers are continuous with respect to 
oo, the above also holds uniformly for oo in suitable small intervals. 

For an implementation we have to discuss how to chose the parameters 
So and Si in order to yield §(cu) C T+ Sj . For the two cases from above (see 
Fig. [4]) we can proceed as follows: 

1. i(0, 00 ) C r+ (Fig. 4(a)): The choice |so| = |si| is sufficient 

Then T 


In partic- 

S0!Sl is a straight line i(so + si) 1 ® 


ular we can chose sq = s 1 . 

Choosing K((sq + si)/|sq + Si|) small enough yields §(cu) C T+ ; 


Sl' 


2. i(—0,0) U i(9, 00 ) C r+ (Fig. 4(b)): We chose arbitrary s 0 ,si, such 


that (3.12) is satisfied. E.g. we can chose so,s~i as in Fig. 4(b)[ Then 
we scale s 0 := s 0 9/((s 0 , sq) and Si := Si9/((s 0 , sq), such that iM fl 
Tso.si = {0, ±i#}. The real valued outgoing wavenumbers are now 
contained in —iF+ S] . To check if the non-real outgoing wavenumbers 
are contained in — vT+ , g SOtSl (in n ) < 1 can be evaluated. If not, define 
Sq := s 0 9/C(s 0 ,ts^ + (1 - t)si) and s\ := (fsh + (1 - t)s 1 )9/((s 0 ,ts^ + 
(1 — t)s 1 ) for t G [0,1]. There hold Sq = s 0 , s? = Si and sj = So, s) = 
So, i.e. t 1 —y r s t s t is a homotopy between r s0;Sl and iM, such that 
r s t s t fl iM = {0, ii#} for all t G [0,1). Choosing t G [0,1] close enough 
to one, S(tu) C t can be ensured. 


The approximation error of the method depends on g SOtSl {\n n ) Nl °^. Hence, 
in general one should try to find s 0 and Si such that g So , Sl ( i« n ) is as small as 
possible for the first wavenumbers n n . 

Note, that the wavenumbers and modes are not used in the numerical 
method. A rough estimate of the location of the wavenumbers is needed in 
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K(lK„) 

(a) w £ (0.8,1.57), sq = Si = —1 + 0.2i 



SR(iK„) 

(b) u £ (1.57,1.78), s 0 = -0.3742-0.48861, 
Si = -0.7752 + 1.0396i 


Figure 4: Lowest outgoing (red) and incoming (blue) wavenumbers (mul¬ 
tiplied with i) for different intervals of frequencies. Separating curves T 
with parameters so (and si) are marked green. Material parameters are 
H = p = E = 1 and v = 0.25. 


order to chose reasonable values for so and si, but the Hardy space infinite 
element method itself is independent of the wavenumbers and the waveguide 
modes. 


4 Tensor product discretization of one waveg¬ 
uide 


For a discretization of (2.1) a standard finite element method for the bounded 
interior domain r2 int can be used. For the waveguides a tensor product 
method appears reasonable. We present here only shortly such a tensor prod¬ 
uct method for one reference waveguide. The extension to several waveguides 
and the coupling with the finite element method for fi int is straightforward. 


We consider again the reference waveguide W := R + x (—R,R). For 
sufficiently fast decaying test functions v we obtain as in the last section a 
variational formulation of (2.4) with the pole condition as radiation condi¬ 
tion: We are looking for u with £u 6 [H~(T S0 ^ S1 ) 0 L 2 (— R, R)] 2 such that 
for all suitable v 


a ext (u,v) — a; 2 6 ex t(u, v) 



W u ) (o)) ’ vds 


(4.1a) 
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with a ext (u ext , v cxt ) := j w (2/ie(u) : e(v) + Adivu divv) d(^,p) and6 ext (u, v) : = 
f w p u • v d(£, rj). In full detail we have 


a ext (u, v) = / ((2 n + A) <%ui <9 ? vi + pd^ux d v vx) d(£, p) + ((2// + A) 0„u 2 S r/ v 2 + /id ? u 2 5 


'TV 


MV 


'TV 


(M,Ui <%v 2 + A^ui <9, ? v 2 ) d(£, rj)+ (pd ? u 2 + A^u 2 <%Vi) d(f, 77 ), 




6ext(u,v) := / p(ui V! + u 2 v 2 ) d(£,r}). 

Jw 


The right hand side of (4.1a) cancels out with the corresponding term of 


the variational formulation of the interior problem, if the test functions are 
chosen to be continuous on the interface T := {0} x (—R,R). 

The whole waveguide is considered as an (in)hnite element with element 
matrix A ext — uj 2 B ext and tensor product basis functions of the form <p* ong <E> 
pti-anS' ]y[ ore precisely, for u we use the ansatz 

jylong ^ytrans 

u(£,p) » Ujviong^r.n.(£,r/) := ^ Ons (0^ rans (^) 

j =1 '=1 


jl 


( 2 ) 

“V 


with aff, aff E C. 

In order to ensure continuity of the discrete solution at the interface T 
the basis functions are the non-vanishing traces of the interior basis 


(€,V) e (0, 00 )x(-R,R) 
(4.2) 


functions E H l (M m t ), i.e. for (p“p| T , • • ■, )I T ^ 0 


'3 

trans 


trans 


T(i) 

vr'~(v) ■= v), v e (-R,R), l = i,...,N 

The surface matrices M trans , D trans , S trans € ^v trans xv trans ( ] e p ne( -] by 

rR 


M t 


trans 


„trans , „trans 


<Pl <Pr, 


dp 


'-R 


N tT 


D 


trans 


l,m= 1 


r R x ^ 

I O trans trans 

' °riVl ‘Pm 
-R 


1 dp 


Si 


trans 


l,m= 1 

(4.3) 

can be computed e.g. by Gauss quadrature. 

For the unbounded ^-direction we use the basis functions and element 
matrices M\ ong , D\ ong and Si on g of the last section (we have droped the su- 


f 

'-R 


perscripts for convenience). The system matrix of (4.1a) is therefore given 


by AL ext - cu 2 B 


2 Dext 


^2Af long ./V trans x 27V long jV trans 


with 


^ext f T A)S'long ® Mtrans T pM\ ong $5 St rans d ^\ong ® -^trans T A -Dlong -D trans 

V d -Diong D trans + A D long $5 D trans (2p + A)Mi ong <8) Strans T P S'long $$ Mfa 

(4.4a) 


trans 


gext ._ I y ±v± long 


P M\ ong $5 Mx rans 0 

0 p M\ ong ® M trans 


(4.4b) 
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The coupling of the matrix A ext — uj 2 B ext with a system matrix A mt — 
l u 2 B mt of the interior problem becomes natural due to <^ ong (0) = 8\j. There¬ 
fore, we have constructed an (in)hnite element for each waveguide which can 
be easily combined with the finite element discretization of the interior prob¬ 
lem. In more details and in 3d this construction is shown in [23j Sec. 4.1], 


5 Resonance problems 


Up to now, we have only considered diffraction problems. We now focus 
on the resonance problem (2.2). Since the modal radiation condition from 


Sec. 2.2 is (more ore less well) defined only for positive frequencies u>, so 


is the resonance problem (2.2). A meaningful extension to complex valued 
frequencies consists in a holomorphic extension of (2.2). If S(a? 


C T + 

so,Si 


holds true for all u in a real interval (a, 6), then the pole condition (2.8) with 
respect to the parameters s 0 , Si is equivalent to the modal radiation condition 
on this interval. Since the pole condition is independent of the frequency, 
it is straightforward to holomorphically extend the resonance problem from 
the interval (a, b ) into the complex plane through the pole condition. 

However, the dependence of this extension on s 0 , Si is not obvious. The 
crucial property to investigate is the Fredholmness of the holomorphic exten¬ 
sion. Due to similar results for acoustic waveguides in [15], we conjecture that 


the essential spectrum, i.e. the set of frequencies u> for which (2.2) equipped 
with the pole condition is not Fredholm, takes the form 




5 (s 0 , Si) := (w 6 C: i/s; n (a;) £ r s0iSl for all n E N} . 


(5.1) 


For complex frequencies to the wavenumbers K n (u) are thereby defined in the 
same way as for real valued frequencies, i.e. as eigenvalues of a quadratic 
eigenvalue problem (see Rem. 2.1). The set <r ess (so,si) can be interpreted as 


the union of N curves, i.e. a ess (s 0 ,Si) = |J neN {u; e C: in n (u) E r s0jS1 }. Due 
to the point symmetry of r s0jSl and of the set of wavenumbers with respect 
to 0, each curve starts at a frequency u, such that there exists a vanishing 
wavenumber k h (uj) = 0. 

To support these conjectures, we give a numerical example in Fig. [5] for 
a waveguide M. x (—1,1) with p = E = 1 and v = 0.25. For a given k E 


—ir s0)S1 the corresponding u> for which wavenumbers K n (oj) < 
are the square roots of a linear eigenvalue problem given in 


-tr 


exist, 
Moreover, 


SO, 51 


-ir 


*0,Sl 


can be explicitly parameterized with — i7 So , Sl due to (3.11). The set 
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Figure 5: Eigenvalues of (2.2) for = M. x (—1,1), p — E — 1 and u = 0.25 


compared with the essential spectrum <r ess (so, Si). 


cr ess (so, si) can therefore be numerically computed with a standard eigenvalue 
solver. For our example we chose 1000 sample points n = — ^7s 0 ,si( r ) £ 
ir so>Sl with equidistant r e (—5,0) to compute cr ess (s 0 , Si). In Fig. [5] we 
compare it with the computed eigenvalues of (2.2), where we used a high 


order hnite element method for h2i nt = (0,5) x (—1,1) with a triangular 
mesh with maximal meshwidth h — 0.1 and polynomial order p = 12. The 
two waveguides (5, oo) x (—1,1) and (—oo, 0) x (—1,1) are discretized using 
Hardy space infinite elements with jV long = 200. For the generalized linear 
eigenvalue problem we have used a shift-and-invert Arnoldi algorithm with 
four different shifts and a Krylov space with dimension 2000. 

It can be seen in Fig. [5] that almost all of the computed eigenvalues A 6 C 
fit very well to cr ess (so, Si). Thus their interpretation as the discretization of 


an essential spectrum is reasonable. The isolated point in Fig. 5(a) has to 


be interpreted as a resonance. We clearly see in Fig. [5] that different pa¬ 
rameter choices so, si lead to different essential spectra, as predicted through 
cr ess (s 0 , si)- We also observe a dependence of the discrete spectrum on the 
parameters. This is because the different parameters correspond to different 
choices of branches of the solution operator. This behavior will be discussed 
in more detail in Sec. 16.31 

Moreover, we observe in Fig. [5] that some curves of the essential spectrum 
are partially located in the quadrant {z G C: ^(z),^ (z) > 0}. Each such 
curve admits in addition to its starting point a supplementary intersection 
uj* with the real axis. This means that there exists a wavenumber n n {u:) 
crossing over —ir s0iSl as cu moves along the real interval ( u * — e,u* + e) 
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(for some e > 0). Since the pole condition selects wavenumbers in — ir+ , 
it can’t be equivalent to the modal radiation condition for both intervals 
(to* — e,co*), (to*, to* + e) simultaneously. A inspection of such cases shows, 
that for intervals / with <T eS s(so,Si) D {x + iy: x G I,y > 0} 7 ^ 0 the pole 
condition is never equivalent to the modal radiation condition. Note, that 
such an interval exists in Fig. |5(a)| indicating the presence of a backward 
propagating mode for ca G (1.58,1.65) (compare with Fig. [3]). Since the 
Hardy space infinite element method only depends on s 0 an d Si, for this 
finding no calculation of wavenumbers or waveguide modes was needed. 

It is worth mentioning that the discrete and essential spectra can also 
be interpreted in terms of the solution operator to ( 2 . 2 ) with respect to the 
pole condition. Resonances correspond to poles and the essential spectrum 
corresponds to branch cuts of the solution operator. 

We conclude, that the computed eigenvalues of (2.2) have to be inter¬ 
preted carefully. They might be part of the discretization of an essential 
spectrum. If this is not the case, they are approximations to the resonances 
of (2.2). These resonances depend on the parameter choice and belong to 
different Riemann sheets. The relevance of resonances for scattering prob¬ 
lems will be investigated further in Sec. 6.2 Note, that the discretization of 
(2.2) with Hardy space infinite elements leads to a generalized linear eigen¬ 
value problem, which can be solved using a standard linear eigenvalue solver. 
This is one of the advantages of the presented method over classical modal 
methods. 


6 Numerical results 

Since rigorous convergence results are not available for Hardy space elements 
in the context of elastic waveguide problems, we report here on numerical 
experiments. The convergence results for acoustic waveguide problems in [T5j 
and of the model problem in [13] indicate a super-algebraic convergence rate 
with respect to the number of unknowns in the infinite direction. The first nu¬ 
merical experiment for a diffraction problem with a known solution supports 
this conjecture. The second numerical experiment underlines the relation be¬ 
tween a scattering problem and the corresponding resonance problem. The 
last experiment illustrates the dependency of the computed resonances on 
the method parameters. 

The computations were all made with the finite element package net- 
gen/NGSolve [26] [27] and the module ngs-waves [23] containing the source 
code for Hardy space methods. For the resonance problems we have used 
a standard shift-and-invert Arnoldi method using MUMPS or PARDISO as 
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Figure 6: Left panel: outgoing (red) and incoming (blue) \K n for c o = 1.66 
. Right panel: relative (Lf 1 (r2; n t)) 2 -error of u — u re f w.r.t. the number of 
unknowns iV long in the infinite direction for different polynomial orders. 

direct solver. 


6.1 Convergence test 

For unperturbed waveguides W = R + x (—i?, R ) the Rayleigh-Lamb modes 
are known to be (semi-)explicit solutions to (2.4) (see e.g. [2]): We define the 
longitudinal wave speed by cl \/(A + 2 p)/p, the transversal wave speed 
by Ct := \/p/p), and for k G C the dispersion relations 


F s (k) := 4:K, 2 a/3 sin(a:R) cos (/3R) + {k 2 — (3 2 ) 2 cos (aR) sin (/3R), (6.1a) 
Fa(k) := 4K 2 cyd cos (aR) sin (/3R) + (k 2 — f3 2 ) 2 sin (aR) cos (/3R), (6.1b) 




(b) error w.r.t. to iV long for p = 10 


Figure 7: relative (7/ 1 (f2 int )) 2 -error of u — u ref for different frequencies u> 
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with a : = a Juj 2 /c\ — k 2 and /3 := \J u 2 jc\ — k 2 . If k is a root of Fs(n) = 0 
or Fa(k) = 0, then 

U S K, v) := e“ 5 w s (,,). := e^w^fa), (£,d) 6 W, 


with 

Sf \_ f ik(k 2 —/3 2 ) 2 sin(/3R) cos(ari)+/32iKasin(aR) cos(/3r]) \ 

^ \V ) • y —a(K 2 —/3 2 ) 2 sin(/3R) sin(ar))+2K 2 asm(aR) sin(/3ri) J ’ 

A( \ _ f iK(K 2 —/3 2 ) 2 cos(/3R)sin(ar))—l32iKacos(oiR)sin(l3r])\ 

^ \V ) y c*(k 2 —/J 2 ) 2 cos(/3_R) cos(ar?)+2K 2 a cos(aR) cos(/3ri) J ’ 


(6.2a) 

(6.2b) 


are the symmetric and anti-symmetric Raylcigh-Lamb modes. 

We chose for R — 1, Young’s modulus E — 1, density p = 1 and Poisson’s 
ratio z/ = 0.25 the reference function 


Uref • — 


S u ?/ll w flU a (T) + Xl u //Il w /IU 2 (T) 


J=1 


J=1 


(6.3) 


with the first five symmetric and first four anti-symmetric Raylcigh-Lamb 
modes. The domain is given by Q — flintUYUbP with h2 int = (0,15) x (—1,1) 
triangulated with maximal mesh size h = 0.25, T = {15} x (—1,1), and 
W = (15, oo) x (—1,1). (2.4) was completed with the Dirichlet boundary 
condition u(0, •) = u re f(0, •) and the pole condition (2.8) for r so , Sl defined 
in Sec. with s 0 = -0.374158 - 0.488609i and s, = -0.775234 + 1.039621. 

First we pick a fixed frequency u j = 1.66, such that there exists an outgo¬ 
ing wavenumber k(uj) < 0 (see Fig. 6(a)) and vary the polynomial order of 
the finite element method for Q int and the number of unknowns N long in the 
infinite direction (see Fig. 6(b)). Super-algebraic convergence in jV long can 
be observed until the error of the finite element discretization of the interior 
problem, which depends on the polynomial order, is reached. 

In Fig. 7(a) we fixed the uniform polynomial degree to p = 5 and varied 
l o G (1.57,1.78). At the frequencies u & 1.6260, u ~ 1.7206 the method fails, 
as already mentioned in Sec. |3.3[ since for these frequencies one outgoing 


wavenumber coincides with an incoming one. Apart from these frequencies, 
the relative error in h2i nt is small. But if we chose to = 1.7205 such that 
there exists one wavenumber n n 0 , the convergence rate w.r.t. Ad ong is 
very poor 


see 


7(b) with p = 10). Techniques to overcome this problem were 


developed for acoustic waveguide problems in [15]. The application to elastic 
waveguide problems is intended for future research. 


6.2 Cavity resonances 

We consider a waveguide Kx (—1,1) with a cavity 9! := (—5, 5) x [1, 6) at¬ 
tached to some wall at the homogeneous Dirichlet-boundary Otto = (—5,5) x 
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{6}. On all other parts of the boundary we assume homogeneous Neumann 
boundary conditions. In Fig. 8(a) the triangulation of the interior domain 


n int = (—10,10) x (—1,1) U fli is given. 

First, we have solved scattering problem (2.1) for different frequencies 
G (1.58,1.65) with an incoming wave consisting of propagating Lamb 


u 


modes and measured the stress in the cavity Qi, which is plotted in Fig.|8(b) 


We have chosen the material parameters as in the last subsection, in order to 
ensure the existence of backward propagating modes in the stated frequency 
interval. 

Second, we solve the resonance problem (2.2) for the same material pa¬ 
rameters and look for resonances near the interval (1.58,1.65). In Fig. [9] the 
computed eigenvalues for two different pairs of s 0 and Si are given. We can 
detect in Fig. [9] four eigenvalues, are separated from the remaining spectrum. 
They are cui = 1.609 — 0.019i, 0 J 2 = 1.625 — 0.003i, cu 3 = 1.655 — 0.003i, and 
oq = 1.66 — O.Olli. Moreover these eigenvalues coincide for both calculations 
obtained with different method parameters. Since the essential spectrum 
depends on the method parameters, we can identify them as resonances. 



Figure 8: scattering of a wave signal by a cavity with Dirichlct boundary 
conditions 


Re to 



Re to 



Figure 9: The computed spectrum for a pair or parameters s 0 , Si indicated as 
red points, and for a slightely perturbed pair s 0 , si indicated as blue squares. 
The right plot is a magnification of one part of the left plot. 
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*(«>„) 

Figure 10: the computed spectrum for s 0 = —0.182086 — 0.2377841, = 

— 1.71492 + 2.29978i and varying a G (5,10 10 ); eigenvalues near to shaded 
domains are discretizations of two curves {w 6 C: i K n (uj) G r s0)Sl } of the 
essential spectrum (see Sec. [5j) with n n being the wavenumbers to one anti¬ 
symmetric Lamb mode and to one symmetric, backward propagating Lamb 
mode. The squares indicate the resonances for a — 11. 


If we compare Fig. 8 (b)| with Fig. [9j we observe that the two peaks in 
Fig. 8(b) fit very well to the two computed resonances ui res from Fig. [9] with 
|7y(u; r es)| most low, i.e. co 2 and 0 J 3 . This is not surprising, since apart from 
the essential spectrum we expect the solution operator to be meromorphic 
with respect to the frequency with resonances being the poles. Hence, in the 
neighborhood of a resonance the scattering problem is almost singular and 
thus very sensitive to external forces. 


6.3 Dependency of resonances with respect to so,si 

In this subsection we want to highlight that the choice of the parameters 
So, Si for solving a resonance problem corresponds to a choice of the Riemann 
sheet on which resonances are sought. We construct an example problem 
parameterized by a parameter a G M + , such that we have some apriori 
knowledge about the resonances and observe the influence of the choice of 
Riemann sheets on the resonances. 

We are solving (2.2) for Q = M x (—1,1) with density p = 4 for x G Hi := 


[—0.5, 0.5] 2 and p = 1 elsewhere, Young’s modulus E = 1 and Poisson’s ratio 
v = 0.2 and B u := er(u) ■ n. Moreover, we added to (2.2) a jump in the 
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Figure 11: real part of first (left panels) and second (right panels) Cartesian 
component of a resonance function for a = 11. 
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normal stress on <9fii 


(cr(u + ) — c(u )) ■ n = au for x G a G 


u mt . v mt _ 


leading to an additional term a J af2 

For the discretization we have used a high order finite element method 
for h2 int = (—2,2) x (—1,1) based on a triangulation with maximal meshsize 
h — 0.1 and a finite element order 10. The waveguides are discretized using 
the Hardy space infinite elements with 150 basis functions in H~(T SOjSl ) and 
so = —0.182086 — 0.2377841, si = —1.71492 + 2.29978i. Since there exists 
a backward propagating mode (see Fig. [2] for uj = 1.615), the parameters 
are chosen such that the pole condition is equivalent to the modal radiation 
condition in the neighborhood of 1.615. 

For a —> oo the additional term leads to two decoupled problems for Hi 
and \ Hi with Dirichlet boundary conditions at <9Qi. Hence, for a —>■ oo 
some of the resonances co a should converge to the square root of the (real) 
Dirichlet-eigenvalues of the problem in the bounded domain f^, which can 
be computed using standard finite element methods. This can be observed 
in Fig. HI The sequences of resonances labeled with 1 and 2 are converging 
for a —» oo to ~ 1.89, which is the square root of a real eigenvalue with 
geometric multiplicity 2. 

In order to distinguish resonances from discretizations of the essential 
spectrum, we have computed as in Sec. [5] additionally the essential spectrum 
of the resonance problem, which depends on so and Si. For a = 11 the real 
part of the Cartesian components of three resonance functions can be seen in 
Fig. 11 We notice that resonances labeled with 1 have symmetric resonance 
functions (as well as these labeled with 3) and the resonances labeled with 2 
have antisymmetric resonance functions. 

In order to illustrate the effect of the essential spectrum on the reso¬ 
nance functions, we have computed first for the sequence 2 of resonances 
Ua in Fig. 12 the wavenumbers K n (iJa ) multiplied with i. The Hardy space 
method computes resonance functions with wavenumbers in — ir + . Since for 
two sequences of symmetric wavenumbers there exists an intersection with 
T and therefore jumps in these wavenumbers, we would expect jumps in the 
resonances 0 Jq \ too. These jumps would be in Fig. [lojexactly at the intersec¬ 
tion points of cja' 1 with the essential spectrum. Since the resonance functions 
corresponding to u are purely antisymmetric, they are orthogonal to the 
symmetric Lamb waves and thus not influenced by the change of symmetric 
wavenumbers. Therefore the resonances labeled with 2 are crossing in Fig. 10 
the curve of the essential spectrum corresponding to a symmetric Lamb mode 
without any perturbation but stop at the curve of the essential spectrum cor- 
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Figure 12: wavenumbers (multiplied with i) for the complex resonances 2 
of Fig. 10 with a G (6,10 10 ) computed with [25], The squares indicate 


the wavenumbers for a = 10 10 of the symmetric Lamb waves, circles of the 


antisymmetric Lamb waves, 
the same as in Fig. 


The parameters for the green curve T 


So,si 


are 


10 


responding to an antisymmetric Lamb mode. The last statement can be also 
seen in Fig. [I2j since for a ~ 6 one antisymmetric wavenumber hits F. 

Fig. [13] shows similar results for the resonances labeled with 1 and 3, where 
the resonance functions are both symmetric. E.g. the resonances (ui 3 ' ) seem 
to vanish after reaching the essential spectrum corresponding to a symmetric 
Lamb mode. The resonances 1 are the counterparts of the resonances 3 on a 
different Riemann sheet: For a = 11 both resonances are present and belong 
to the same Riemann sheet, which is determine d by the chosen parameters So 


and Si. The modal radiation condition of Sec. 2.2 for the real parts 

/o\ 

and } of the resonances support a backward propagating mode (for 

the real parts of the res onanc es see Fig. 10, the corresponding dispersion 
curves are given in Fig. 2(a)). is located above and below the 

essential spectrum arising from a branch cut of the holomorphic extension of 

/o\ 

this backward propagating mode. For a ~ 11.7 the resonances co y a hit the 
branch cut and vanish with increasing a, since the Riemann sheet of their 
existence is hidden in t his p art of the complex plane. Hence, for larger values 


of a we can see in Fig. 10 only the resonances c Ua^ 
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3 



(a) a € (5,10 10 ). The squares indicate the 
wavenumbers for a = 10 10 of the symmet¬ 
ric Lamb waves, circles of the antisymmetric 
Lamb waves 



(b) a e (9.4,11.7) 


Figure 13: wavenumbers (multiplied with i) for the complex resonances 1 
(left) and 3 (right) of Fig. [lOjw.r.t. a computed with [25] . The parameters 
for the green curve r so S1 are the same as in Fig. 10 


7 Conclusion 


In this paper we addressed the topic of backward propagating modes in time- 
harmonic two-dimensional elastic waveguides. The so called pole condition 
allows to reformulate the radiation condition even in the presence of back¬ 
ward waves without using the waveguide modes or wavenumbers. For the 
problems under consideration a special class of pole conditions depending on 
two complex parameters s o and Si is sufficient. Detailed explanations on the 
choice of these parameters are given in Section [373 The Hardy space infinite 
element method developed in [I3J relies on this pole condition and can be 
easily applied to waveguide structures using tensor product elements. 

An outstanding property of the presented scheme is that it is independent 
of the waveguide modes and wavenumbers which simplifies the implementa¬ 
tion. Moreover, a discetization of a resonance problem leads to a generalized 
linear matrix eigenvalue problem, which is much easier to handle than a non¬ 
linear problem. However, the interpretation of the occurring spectral objects 
for resonance problems is rather involved and was addressed in Section [5] and 
Section 16.31 


















The method shows in numerical examples super-algebraic convergence 
for diffraction problems. For resonance problems the numerical results also 
coincide with the theoretical considerations. In particular, a strong relation 
between the resonances and the behavior of the solutions to diffraction prob¬ 
lems with frequencies in the neighborhood of a resonance frequency can be 
seen. Hence, resonance and diffraction problems in two dimensional waveg¬ 
uides can be reliably solved using the Hardy space infinite element method. 
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